home *** CD-ROM | disk | FTP | other *** search
/ Enter 2004 January / enter-2004-01.iso / files / maxima-5.9.0.exe / {app} / share / maxima / 5.9.0 / src / numerical / slatec / zseri.lisp < prev    next >
Encoding:
Text File  |  2003-02-09  |  9.3 KB  |  242 lines
  1. ;;; Compiled by f2cl version 2.0 beta 2002-05-06
  2. ;;; 
  3. ;;; Options: ((:prune-labels nil) (:auto-save t) (:relaxed-array-decls t)
  4. ;;;           (:coerce-assigns :as-needed) (:array-type ':simple-array)
  5. ;;;           (:array-slicing nil) (:declare-common nil)
  6. ;;;           (:float-format double-float))
  7.  
  8. (in-package "SLATEC")
  9.  
  10.  
  11. (let ((zeror 0.0) (zeroi 0.0) (coner 1.0) (conei 0.0))
  12.   (declare (type double-float conei coner zeroi zeror))
  13.   (defun zseri (zr zi fnu kode n yr yi nz tol elim alim)
  14.     (declare (type (simple-array double-float (*)) yr yi)
  15.              (type f2cl-lib:integer4 kode n nz)
  16.              (type double-float zr zi fnu tol elim alim))
  17.     (prog ((wr (make-array 2 :element-type 'double-float))
  18.            (wi (make-array 2 :element-type 'double-float)) (i 0) (ib 0)
  19.            (idum 0) (iflag 0) (il 0) (k 0) (l 0) (m 0) (nn 0) (nw 0) (aa 0.0)
  20.            (acz 0.0) (ak 0.0) (ak1i 0.0) (ak1r 0.0) (arm 0.0) (ascle 0.0)
  21.            (atol 0.0) (az 0.0) (cki 0.0) (ckr 0.0) (coefi 0.0) (coefr 0.0)
  22.            (crscr 0.0) (czi 0.0) (czr 0.0) (dfnu 0.0) (fnup 0.0) (hzi 0.0)
  23.            (hzr 0.0) (raz 0.0) (rs 0.0) (rtr1 0.0) (rzi 0.0) (rzr 0.0) (s 0.0)
  24.            (ss 0.0) (sti 0.0) (str 0.0) (s1i 0.0) (s1r 0.0) (s2i 0.0)
  25.            (s2r 0.0))
  26.       (declare (type (simple-array double-float (2)) wr wi)
  27.                (type double-float s2r s2i s1r s1i str sti ss s rzr rzi rtr1 rs
  28.                 raz hzr hzi fnup dfnu czr czi crscr coefr coefi ckr cki az atol
  29.                 ascle arm ak1r ak1i ak acz aa)
  30.                (type f2cl-lib:integer4 nw nn m l k il iflag idum ib i))
  31.       (setf nz 0)
  32.       (setf az (zabs zr zi))
  33.       (if (= az 0.0) (go label160))
  34.       (setf arm (* 1000.0 (f2cl-lib:d1mach 1)))
  35.       (setf rtr1 (f2cl-lib:fsqrt arm))
  36.       (setf crscr 1.0)
  37.       (setf iflag 0)
  38.       (if (< az arm) (go label150))
  39.       (setf hzr (* 0.5 zr))
  40.       (setf hzi (* 0.5 zi))
  41.       (setf czr zeror)
  42.       (setf czi zeroi)
  43.       (if (<= az rtr1) (go label10))
  44.       (multiple-value-bind
  45.           (var-0 var-1 var-2 var-3 var-4 var-5)
  46.           (zmlt hzr hzi hzr hzi czr czi)
  47.         (declare (ignore var-0 var-1 var-2 var-3))
  48.         (setf czr var-4)
  49.         (setf czi var-5))
  50.      label10
  51.       (setf acz (zabs czr czi))
  52.       (setf nn n)
  53.       (multiple-value-bind
  54.           (var-0 var-1 var-2 var-3 var-4)
  55.           (zlog hzr hzi ckr cki idum)
  56.         (declare (ignore var-0 var-1))
  57.         (setf ckr var-2)
  58.         (setf cki var-3)
  59.         (setf idum var-4))
  60.      label20
  61.       (setf dfnu (+ fnu (f2cl-lib:int-sub nn 1)))
  62.       (setf fnup (+ dfnu 1.0))
  63.       (setf ak1r (* ckr dfnu))
  64.       (setf ak1i (* cki dfnu))
  65.       (setf ak
  66.               (multiple-value-bind
  67.                   (ret-val var-0 var-1)
  68.                   (dgamln fnup idum)
  69.                 (declare (ignore var-0))
  70.                 (setf idum var-1)
  71.                 ret-val))
  72.       (setf ak1r (- ak1r ak))
  73.       (if (= kode 2) (setf ak1r (- ak1r zr)))
  74.       (if (> ak1r (- elim)) (go label40))
  75.      label30
  76.       (setf nz (f2cl-lib:int-add nz 1))
  77.       (f2cl-lib:fset (f2cl-lib:fref yr (nn) ((1 n))) zeror)
  78.       (f2cl-lib:fset (f2cl-lib:fref yi (nn) ((1 n))) zeroi)
  79.       (if (> acz dfnu) (go label190))
  80.       (setf nn (f2cl-lib:int-sub nn 1))
  81.       (if (= nn 0) (go end_label))
  82.       (go label20)
  83.      label40
  84.       (if (> ak1r (- alim)) (go label50))
  85.       (setf iflag 1)
  86.       (setf ss (/ 1.0 tol))
  87.       (setf crscr tol)
  88.       (setf ascle (* arm ss))
  89.      label50
  90.       (setf aa (exp ak1r))
  91.       (if (= iflag 1) (setf aa (* aa ss)))
  92.       (setf coefr (* aa (cos ak1i)))
  93.       (setf coefi (* aa (sin ak1i)))
  94.       (setf atol (/ (* tol acz) fnup))
  95.       (setf il (min (the f2cl-lib:integer4 2) (the f2cl-lib:integer4 nn)))
  96.       (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
  97.                     ((> i il) nil)
  98.         (tagbody
  99.           (setf dfnu (+ fnu (f2cl-lib:int-sub nn i)))
  100.           (setf fnup (+ dfnu 1.0))
  101.           (setf s1r coner)
  102.           (setf s1i conei)
  103.           (if (< acz (* tol fnup)) (go label70))
  104.           (setf ak1r coner)
  105.           (setf ak1i conei)
  106.           (setf ak (+ fnup 2.0))
  107.           (setf s fnup)
  108.           (setf aa 2.0)
  109.          label60
  110.           (setf rs (/ 1.0 s))
  111.           (setf str (- (* ak1r czr) (* ak1i czi)))
  112.           (setf sti (+ (* ak1r czi) (* ak1i czr)))
  113.           (setf ak1r (* str rs))
  114.           (setf ak1i (* sti rs))
  115.           (setf s1r (+ s1r ak1r))
  116.           (setf s1i (+ s1i ak1i))
  117.           (setf s (+ s ak))
  118.           (setf ak (+ ak 2.0))
  119.           (setf aa (* aa acz rs))
  120.           (if (> aa atol) (go label60))
  121.          label70
  122.           (setf s2r (- (* s1r coefr) (* s1i coefi)))
  123.           (setf s2i (+ (* s1r coefi) (* s1i coefr)))
  124.           (f2cl-lib:fset (f2cl-lib:fref wr (i) ((1 2))) s2r)
  125.           (f2cl-lib:fset (f2cl-lib:fref wi (i) ((1 2))) s2i)
  126.           (if (= iflag 0) (go label80))
  127.           (multiple-value-bind
  128.               (var-0 var-1 var-2 var-3 var-4)
  129.               (zuchk s2r s2i nw ascle tol)
  130.             (declare (ignore var-0 var-1 var-3 var-4))
  131.             (setf nw var-2))
  132.           (if (/= nw 0) (go label30))
  133.          label80
  134.           (setf m (f2cl-lib:int-add (f2cl-lib:int-sub nn i) 1))
  135.           (f2cl-lib:fset (f2cl-lib:fref yr (m) ((1 n))) (* s2r crscr))
  136.           (f2cl-lib:fset (f2cl-lib:fref yi (m) ((1 n))) (* s2i crscr))
  137.           (if (= i il) (go label90))
  138.           (multiple-value-bind
  139.               (var-0 var-1 var-2 var-3 var-4 var-5)
  140.               (zdiv coefr coefi hzr hzi str sti)
  141.             (declare (ignore var-0 var-1 var-2 var-3))
  142.             (setf str var-4)
  143.             (setf sti var-5))
  144.           (setf coefr (* str dfnu))
  145.           (setf coefi (* sti dfnu))
  146.          label90))
  147.       (if (<= nn 2) (go end_label))
  148.       (setf k (f2cl-lib:int-sub nn 2))
  149.       (setf ak (coerce (the f2cl-lib:integer4 k) 'double-float))
  150.       (setf raz (/ 1.0 az))
  151.       (setf str (* zr raz))
  152.       (setf sti (* (- zi) raz))
  153.       (setf rzr (* (+ str str) raz))
  154.       (setf rzi (* (+ sti sti) raz))
  155.       (if (= iflag 1) (go label120))
  156.       (setf ib 3)
  157.      label100
  158.       (f2cl-lib:fdo (i ib (f2cl-lib:int-add i 1))
  159.                     ((> i nn) nil)
  160.         (tagbody
  161.           (f2cl-lib:fset (f2cl-lib:fref yr (k) ((1 n)))
  162.                          (+
  163.                           (* (+ ak fnu)
  164.                              (-
  165.                               (* rzr
  166.                                  (f2cl-lib:fref yr
  167.                                                 ((f2cl-lib:int-add k 1))
  168.                                                 ((1 n))))
  169.                               (* rzi
  170.                                  (f2cl-lib:fref yi
  171.                                                 ((f2cl-lib:int-add k 1))
  172.                                                 ((1 n))))))
  173.                           (f2cl-lib:fref yr ((f2cl-lib:int-add k 2)) ((1 n)))))
  174.           (f2cl-lib:fset (f2cl-lib:fref yi (k) ((1 n)))
  175.                          (+
  176.                           (* (+ ak fnu)
  177.                              (+
  178.                               (* rzr
  179.                                  (f2cl-lib:fref yi
  180.                                                 ((f2cl-lib:int-add k 1))
  181.                                                 ((1 n))))
  182.                               (* rzi
  183.                                  (f2cl-lib:fref yr
  184.                                                 ((f2cl-lib:int-add k 1))
  185.                                                 ((1 n))))))
  186.                           (f2cl-lib:fref yi ((f2cl-lib:int-add k 2)) ((1 n)))))
  187.           (setf ak (- ak 1.0))
  188.           (setf k (f2cl-lib:int-sub k 1))
  189.          label110))
  190.       (go end_label)
  191.      label120
  192.       (setf s1r (f2cl-lib:fref wr (1) ((1 2))))
  193.       (setf s1i (f2cl-lib:fref wi (1) ((1 2))))
  194.       (setf s2r (f2cl-lib:fref wr (2) ((1 2))))
  195.       (setf s2i (f2cl-lib:fref wi (2) ((1 2))))
  196.       (f2cl-lib:fdo (l 3 (f2cl-lib:int-add l 1))
  197.                     ((> l nn) nil)
  198.         (tagbody
  199.           (setf ckr s2r)
  200.           (setf cki s2i)
  201.           (setf s2r (+ s1r (* (+ ak fnu) (- (* rzr ckr) (* rzi cki)))))
  202.           (setf s2i (+ s1i (* (+ ak fnu) (+ (* rzr cki) (* rzi ckr)))))
  203.           (setf s1r ckr)
  204.           (setf s1i cki)
  205.           (setf ckr (* s2r crscr))
  206.           (setf cki (* s2i crscr))
  207.           (f2cl-lib:fset (f2cl-lib:fref yr (k) ((1 n))) ckr)
  208.           (f2cl-lib:fset (f2cl-lib:fref yi (k) ((1 n))) cki)
  209.           (setf ak (- ak 1.0))
  210.           (setf k (f2cl-lib:int-sub k 1))
  211.           (if (> (zabs ckr cki) ascle) (go label140))
  212.          label130))
  213.       (go end_label)
  214.      label140
  215.       (setf ib (f2cl-lib:int-add l 1))
  216.       (if (> ib nn) (go end_label))
  217.       (go label100)
  218.      label150
  219.       (setf nz n)
  220.       (if (= fnu 0.0) (setf nz (f2cl-lib:int-sub nz 1)))
  221.      label160
  222.       (f2cl-lib:fset (f2cl-lib:fref yr (1) ((1 n))) zeror)
  223.       (f2cl-lib:fset (f2cl-lib:fref yi (1) ((1 n))) zeroi)
  224.       (if (/= fnu 0.0) (go label170))
  225.       (f2cl-lib:fset (f2cl-lib:fref yr (1) ((1 n))) coner)
  226.       (f2cl-lib:fset (f2cl-lib:fref yi (1) ((1 n))) conei)
  227.      label170
  228.       (if (= n 1) (go end_label))
  229.       (f2cl-lib:fdo (i 2 (f2cl-lib:int-add i 1))
  230.                     ((> i n) nil)
  231.         (tagbody
  232.           (f2cl-lib:fset (f2cl-lib:fref yr (i) ((1 n))) zeror)
  233.           (f2cl-lib:fset (f2cl-lib:fref yi (i) ((1 n))) zeroi)
  234.          label180))
  235.       (go end_label)
  236.      label190
  237.       (setf nz (f2cl-lib:int-sub nz))
  238.       (go end_label)
  239.      end_label
  240.       (return (values nil nil nil nil nil nil nil nz nil nil nil)))))
  241.  
  242.